Learning the nonlinear interactions from particle trajectories 
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Nonlinear interaction of membrane proteins with cytoskeleton and membrane leads to non- 
Gaussian structure of their displacement probability distribution. We propose a novel statistical 
analysis technique for learning the characteristics of the nonlinear potential from the cumulants of 
the displacement distribution. The efficiency of the approach is demonstrated on the analysis of 
kurtosis of the displacement distribution of the particle traveling on a membrane in a cage-type 
potential. Results of numerical simulations are supported by analytical predictions. We show that 
the approach allows robust identification of the potential for the much lower temporal resolution 
compare with the mean square displacement analysis. 



Rapid progress in video-capturing, sub-diffractive mi- 
croscopes and fluorescence technologies has transformed 
a single particle tracking (SPT) technology in a pow- 
erful tool for studying the properties of biological envi- 
ronments and complex fluids [1]. In typical experiment 
some specific type of biomolecule, i.e. protein or lipid 
is labeled by fluorophore or nanoparticle, and its mo- 
tion is tracked with a camera through in subdiffractive 
resolution [3] . The abundance of data available from the 
SPT experiments has risen demand in data analysis tech- 
niques that would help scientists to characterize the in- 
teraction of particle with the environment based on the 
statistical properties of particle trajectories. Most of the 
currently used approaches rely on the analysis of second 
order moments, like mean square displacement (MSD). 
The main objective of this Letter is to demonstrate the 
potential of other statistical properties, that go beyond 
Gaussian approximation and second-order correlations. 
In a practically interesting example of protein moving on 
the membrane we show that many characteristics of the 
particle-membrane interactions that can not be recovered 
from the analysis of MSD reveal themselves as distinct 
statistical signatures in the kurtosis of the particle dis- 
placement distribution. 

The motion of proteins and lipids within biological 
membranes plays important role in many biological pro- 
cesses. Previous assumption that biological membrane 
can be considered as two-dimensional fluid with freely 
diffusing lipids and proteins [2, 3 J are now significantly al- 
tered by the experimental observations that membranes 
are highly heterogeneous [5J [Hj- Models of lipid rafts, 
pickets and fences, protein-protein complexes and protein 
islands were suggested [7HT0]. According to these mod- 
els the motion of the location and diffusion of membrane 
proteins are significantly influenced by the domains of 
different lipid or protein compositions (lipid rafts, pro- 
tein islands) as well as by the interaction with the cy- 
toskeleton and anchored transmembrane proteins (form 
fences and pickets, respectively). E.g., the compartmen- 
talization of the plasma membrane is perhaps the best 
explained by the fences and pickets [H] . Inside each com- 



partment proteins (lipids) experience fast diffusion (at 
time scales < 0.01s [S|) which agrees well with the diffu- 
sion in the artificial membranes (which do not have cy- 
toskeleton). A hopping between different compartments 
(jump over fences) occurs at larger time scale t% ~ 0.01s. 
(see e.g. Table 2 in Ref. [8] for the specific values of tu 
for several types of cells). 

Most SPT studies have been relying on the standard 
video rate (~ 30 frames/sec) (T] which does not allow 
a detailed resolution of the fast diffusion inside com- 
partments because inter-compartment hoping rate I/77, 
exceeds the video rate. The exception is the work of 
Kusumi group [5] which uses 25/is temporal resolution. 
The analysis of SPT trajectories in the vast majority of 
previous work has been based on the analysis of MSD pQ . 
It was demonstrated that SPT with the standard video 
rate is not sufficient to recover any details about fast dif- 
fusion inside compartments as well as any information 
about compartments |S]. MSD uses only a small part 
of information about properties of particle trajectories. 
The only exception when MSD is optimal corresponds 
to the pure random walk of the particle when proba- 
bility distribution of particle displacement is Gaussian. 
However, any inhomogeneity on plasma membrane (rep- 
resented e.g. by the inhomogeneous potential) results 
in the non-Gaussianity of that probability distribution 
which makes MSD non-optimal to recover the properties 
of the system from the SPT trajectories. 

In parallel to biology SPT based approaches were de- 
veloped in microfluidics under the name of microrheology 
|12) . Tracking of the Brownian motion of individual par- 
ticles immersed in viscoelastic fluid, allows reconstruc- 
tion of viscoelastic modulus from the Laplace transform 
of particle MSD. To our knowledge, all common varia- 
tions of microrheology (including two-particle and active 
microrheology) are based on the analysis of second-order 
correlation functions, and assume linear response of the 
viscoelastic fluid. Although microrheological settings are 
not described in this Letter, the methods discussed below 
can be naturally applied there. 

In this Letter we propose to recover major features of 
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the potential from SPT trajectories using the kurtosis as 
the measure of non-Gaussianity. We demonstrate by the 
combination of numerical and analytical methods that 
the rate ~ 100 frames/sec might be sufficient for that 
purpose far superior to the performance of MSD-based 
methods. We are not aimed to fully recover the spatial 
distribution of the potential because such type of inverse 
problem is ill-posed pTJ and requires very large statistics 
of SPT trajectories not available in experiment. 

The starting point of our work is the following obser- 
vation. Whenever a particle experiences nonlinear inter- 
actions, for example with the cytoskeleton, the proba- 
bility distribution of the particle displacement becomes 
non-Gaussian. Therefore, analysis of the deviations from 
Gaussian distribution can reveal new information about 
the nonlinear interactions. There is a infinite num- 
ber of characteristics that measure the degree of non- 
Gaussianity, because the nonlinear interaction can take 
infinite number of forms. In this Letter we consider only 
one of the characteristics, that can be accurately esti- 
mated with limited amount of time-series data. Kurtosis 
of the particle displacement distribution, defined in the 
following way. 
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Here Ar a (t) = r a (t) — r a (0) is the a component of the 
particle displacement, and the angular brackets denote 
averaging over the ensemble of particle trajectories in 
the same nonlinear environment. Whenever the dynam- 
ics of the particle is described by linear equations, the 
resulting distribution will be Gaussian, and the kurtosis 
will be equal to zero for any time interval t. It is there- 
fore very different from the extensively studied anoma- 
lous diffusion described by MSD (Ar a (t) 2 ) that can be 
observed even if the dynamics is linear and the distri- 
bution is Gaussian. Kurtosis of the displacement distri- 
bution originates from the nonlinear interactions of the 
particle with the environment, and as will be shown be- 
low, incorporates a lot of information about the structure 
of these interactions. 

For modeling purposes we consider the motion of 
Brownian particle (biomolecule) on a membrane with co- 
ordinate r = (x,y) in a potential U(r) governed by the 
overdamped Langevin equation 



dt 
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where D is the diffusion coefficient, 7 is the friction co- 
efficient and is the random noise with the Gaussian 
statistics, zero mean and unit variance so that \/2Z?£(i) 
represents the effect of the thermal noise. 

The equation ^ is simulated using the Monte Carlo 
algorithm as a two-dimensional (2D) random walk on 
an elementary triangular lattice composed of equilateral 
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FIG. 1: Simulation of the trajectory of a diffusing protein 
in triangular compartments of size 600 nm during 10 seconds 
shown with a time resolution 10 ms. The dashed lines indicate 
the borders between compartments. 



triangles with side length a = 1 nm (or a = 0.25 nm in 
the case of the inset in Figure [2J . 

The membrane compartments are typically modeled as 
equilateral triangles filling entire 2D plane with sides of 
length L = 300 nm, 150 nm and 600 nm. The potential 
energy on the elementary lattice is labeled as Ui = U(ri), 
where i is 2D index on that lattice. Figure [l] shows the 
typical simulation of total time 10 s. In each simulation 
step, a random closest neighbor j of site i that is occu- 
pied by the random walker (diffusing protein) is chosen. 
The move to site j from i is accepted with probability 
p = min(l, exp [(Ui — Uj)/T]). We set the temperature 
T = 1 which implies by the Einstein relation D = T/7 
that 7 = 1/D in Q . The potential of barriers between 
compartments is defined as Ui(l) — H exp (— d 2 t /a 2 ), 
where Ui(l) is the the contribution to the total poten- 
tial Ui at site i from Ith barrier and dij is the distance 
from site i to Ith barrier. Ui is given by the sum of con- 
tributions from all barriers: Ui = J2i Ui(l). Also H is 
the height of the barrier and the width of the barrier is 
2cr. By default the barrier parameter a is set to 5 nm 
and the height is set to H = 7. 

We set the diffusion coefficient to D = 2.5/ito 2 s~ 1 on 
the lattice, which means that one iteration step of the 
MC simulation corresponds to a time period r = 0.1/xs 
for lattice size a = lnm and r = 0.00625/xs for a — 
0.25nm according to the relation D = a 2 /(4r). Effect 
of discretization is negligible at time scales 3> t which 
motivates our choice of the numerical values of a. 

In absence of the potential U the particle experiences a 
Brownian motion with (Ax 2 +Ay 2 ) = 4Dt and K(t) = 0. 
However, if U is nonzero, the kurtosis of the particle mo- 
tion acquires a non-trivial shape as shown in Figure [2] for 
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FIG. 2: The kurtosis of the displacements along x for 3 differ- 
ent compartment sizes calculated with the resolution 0.1ms. 
The inset shows kurtosis scaled as 1/L for the resolution 1/xs. 



three different compartment sizes. The typical length of 
simulation was 10 9 steps with every 10 3 -th step recorded, 
thus corresponding to an experiment of total length 100s 
with a time resolution (i.e. inverse frame rate) 0.1 ms. 
Note that we always use the elementary time step t to 
generate particle trajectories. Experimental observations 
have much lower temporal resolution than t and to imi- 
tate such resolution we use only a small fraction of sim- 
ulation points to calculate the kurtosis for each chosen 
resolution. Figure [2] shows that the kurtosis is charac- 
terized by two peaks separated by a local minima. Note 
that the position of the minima scales as L 2 . The inset of 
Figure [2] shows a simulation with lattice size a = 0.25nm, 
the duration 10s and the resolution 1/is, which we used 
to test the analytical predictions for kurtosis behavior for 
low times as discussed below. The kurtosis in that inset 
is divided by L. The kurtosis scales as 1/L for low times 
t and grows as t 1 / 2 , as expected from our analysis. 

The kurtosis of simulations with three different tem- 
poral resolutions 0.1ms, 1ms and 10ms for the compart- 
ment size L = 300nm are shown in Figure [J Already 
for the resolution ~ 10ms it is possible to see the charac- 
teristic features of the kurtosis inferring compartments' 
structure by comparing with the pure diffusion case. The 
kurtosis curves are averaged over 5 different simulation 
runs, each of duration 100s for the resolution 0.1ms and 
500s for the resolutions 1ms and 5ms. We also show 
the kurtosis for pure diffusion simulations, which should 
be equal to for the infinite trajectory. The error bars 
show standard deviation for the five simulations with the 
resolution 0.1ms. The inset compares MSD of pure diffu- 
sion with no barriers and two different simulations with 
the resolutions 0.1ms and 10ms for L — 300nm. While 
the MSD plot with the resolution 0.1ms shows transi- 



FIG. 3: The kurtosis along x for 3 different resolutions with 
L — 300 nm. The kurtosis for pure diffusion (no barriers) is 
also shown with the resolution 0.1ms. Error bars correspond 
to 0.1ms resolution. The inset shows MSD for diffusion with 
no compartments (short dashed line) and for simulations with 
the resolutions 0.1ms (long dashed line) and 10ms (solid line). 
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FIG. 4: The kurtosis along x for different fractions of ran- 
domly removed barriers for simulation of 100s duration. 



tion from the fast diffusion regime inside the compart- 
ment with D = 2.5/im 2 s _1 to the slow hopping diffusion 
regime with D = 0.1/im 2 s _1 , the MSD plot of 10ms res- 
olution is able to capture only the slow diffusion. 

Figure [4] shows that the kurtosis is quite robust to the 
defects of the triangular compartment lattice. We also 
found similar robustness when we randomly distorted 
equilateral triangular lattice by 10% in angles compare 
with 7r/3 angles as well as when we used rectangular lat- 
tice instead of triangular one. 

Observed structure of the kurtosis can be explained 
theoretically in a limiting case when the characteristic 
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time Td ~ L 2 /D associated with the diffusion inside the 
compartment is much smaller than the hopping time that 
estimates the interval between consecutive hopping over 
the barrier: t/j ~ exp(AU/T), where AU is the height 
of the barrier. 

The initial rise of the kurtosis function K{t) in the 
region t <C Td is related to short trajectories that were 
reflected from the boundary. In order to understand this 
rise qualitatively one can analyze a simplified problem of 
one-dimensional diffusion in the neighborhood of reflect- 
ing wall. Assuming that the random walk takes place in 
the x direction, the Green function which is defined as 
a probability density of the final position x — X(t) as- 
suming that the particle was at position x$ at t = is 
given by G(x; xq, t) = Gq(x — xq; t) + G(x + xq; t), where 
G (x;t) = (2nDt)- 1 / 2 exp [-x 2 /2Dt] . In the above we 
assumed that the reflecting boundary is located at x = 0. 
The kurtosis can be found by calculating the expres- 
sions for moments C n = J °° dx(x — xo) n G(x; xq, t) with 
n = 2, 4. The moments have to be then averaged over the 
value of xq, that we assume to be uniformly distributed 
in the region < x < L with L ^> y/Dt being a char- 
acteristic compartment size. Straightforward integration 
yields the following expressions for the moments in the 
leading order over y/Dt/L: C 2 = Dt/2 - (Dt) 3 / 2 /3y/icL 
and C 4 = 3(Dt) 2 /4 - i(Dt) 5 / 2 /5y/^L 2 . Note that the 
correction to C4 is smaller in comparison to the cor- 
rection to Cf which explains the initial rise of kurtosis 
K = A(Dt) 1 / 2 /15y/wL. The numerical factors in this 
expression are not universal and may depend on the ac- 
tual form of the compartment. However, the scaling laws 
K cx t 1 ' 2 and K oc L 1 are universal, they are seen in 
the inset of Figure [2] and can be checked experimentally. 

For intermediate timescales -C t -C 7), the parti- 
cle has enough time to diffuse around its compartment, 
however the events of passing through the barrier are still 
rare. In this regime one can calculate the value of kurtosis 
by assuming that the Green function G(; ro, t) — Poo(r), 
where P^ (r) is the uniform distribution with support in- 
side the compartment. This assumption implies that the 
particle had enough time to explore the whole compart- 
ment, and moreover, that the width of the barriers is 
negligible in comparison to the compartment size. If the 
later approximation is not justified, the Green function 
has to be replaced by equilibrium Boltzmann-type dis- 
tribution G(r;ro,t) — exp[—U(r)/T]/Z. As long as the 
initial particle position is also uniformly distributed over 
the compartment we obtain the following expressions for 
the moments of particle jumps: C n — f c dr J c ro(r — 
r ) n Poo(r)Poo(ro) that yields K(t) = -1/10 for triangu- 
lar compartments with very thin barriers. Note, that the 
value of kurtosis on these timescales is sensitive to the 
actual shape of compartments and to the width of the 
barrier potentials and can therefore be used for getting 
experimental insight on these membrane properties. 

The late time asymptote t > Th, that is responsible 



for the second peak of the observed kurtosis is deter- 
mined by trajectories that had a finite number of hops 
between the compartments. The non-gaussianity of the 
jump distance distribution is related to the discrete na- 
ture of barrier hopping events. When the typical number 
of hopping events is small, say 2 — 3, the fluctuations of 
the total distance traveled by a particle are stronger than 
Gaussian, and that explains the rise of the kurtosis at 
t m Th- As the number of hops becomes very large for 
t Th, the distribution becomes Gaussian again, due to 
central limit theorem. The kurtosis decays back to zero. 
Analytical results for the kurtosis behavior in this region 
will be reported elsewhere. 

To conclude, we have proposed and analyzed a novel 
particle-tracking approach for identification of nonlinear 
interactions. Unlike common MSD techniques, our ap- 
proach is based on the analysis of non-Gaussian charac- 
teristics of the particle dynamics, specifically the kurtosis 
of the displacement probability distribution. The func- 
tional dependence of the kurtosis on the measurement 
time carries a lot of information about the nonlinear in- 
teractions that contribute to the particle motion. E.g., if 
a particle is placed in a cage-type potential induced by cy- 
toskeleton or transmembrane proteins, the resulting kur- 
tosis of the displacement is a non-monotonous function 
with three distinct regions characterized by the change 
of the sign of the kurtosis slope. Specific structure of the 
kurtosis function depends on the characteristics of the 
potential: shape and size of the individual cells, heights 
and widths of the barriers but the change of sign feature 
is quite robust to the specifics of the potential. 

A number of other techniques have been proposed for 
identification of the potentials on biological membranes 
based on the analysis of individual trajectories. Most 
comprehensive approach is to solve the inverse problem 
of reconstruction of U (r) from the trajectories. However, 
such problem is generally ill-posed pj] and needs very 
large statistics of trajectories. E.g., Ref.fTB"] suggested 
to infer forces acting on the biomolecule |14| requiring 
the multiple particle visits of each spatial location which 
is difficult to achieve experimentally. In addition, the 
potential U in living cells can slowly change with time. 
This fact may limit the application of inverse problem 
approaches, that attempt to reconstruct the specific po- 
tential. Another approach [S1[TS] focuses on identification 
of the potential barriers from MSD-based analysis. That 
approach is successful but requires very high temporal 
resolution of trajectories. One more method is based on 
the measurement of autocorrelation of SPT trajectories 
and recovering of the probability distribution of parti- 
cle jumps [16 which indirectly displays the information 
about the inhomogeneity of the plasma membrane. In 
contrary, the technique analyzed in this Letter has more 
modest goals of learning the characteristic scales associ- 
ated with the potential. It does not rely on long observa- 
tions of individual particles, and can be based on aggre- 
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gation of time-series from an ensemble of measurements 
of different particles in the same class of membranes. 

Work of P.L. was supported by NSF grant DMS 
0719895. 
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